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ANL-CANDID, 


A Two-dimensional, Diffusion-theory Code 
Based on CANDIDZD 


by 


G. K. Leaf, A. S. Kennedy, 
and G. C. Jensen 


ABSTRACT 


ANL-CANDID is anextension and modification of the 
two-dimensional diffusion-theory code CANDID2D developed 
by Computer Applications Incorporated’ foruse on the Con- 
trol Data Corporation 3600 computer. The code, as delivered, 
was capable of performing a reactivity calculation as well 
as geometry, composition, and buckling searches. These cal- 
culations were limited to rz geometry. In addition, the code 
could perform an up-scattering calculation. The present 
code is the result of extending and modifying the C.A.I. code. 
The extensions include the capability of performing calcu- 
lations in xy and r@ geometry, including the full periodic 
case. Adjoint and source calculations were added with an 
acceleration procedure for the source calculation. In addi- 
tion, an a calculation was added. The use of two-term 
Chebyshev extrapolation was retained; however, the strategy 
for employing it was changed. In addition, the strategy em- 
ployed in the search procedure was completely changed. 
Mesh refinement and restart capability have been provided. . 


I. FINITE-DIFFERENCE EQUATIONS 


A. Derivation of Equations for r@ Geometry 


We shall begin our discussion by deriving the finite-difference 
equations for r@ geometry. The derivation is based on the usual first- 
order approximation and is consistent with the derivation for rz geometry.’ 
In the absence of external sources, the system of multigroup diffusion 
equations can be written in the form | 


-div(D,V@,) + (08 +D oP )be= > o8'8¢., 
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with external boundary conditions having the form 


A,Vobg “nt Boog = Cy, (152) 


where B denotes the transverse buckling, and n is the unit normal for the 
exterior surface of the reactor. The reactor configuration is assumed to 
be composed of regions such that the macroscopic cross sections are con- 
stant within each region. 


The object of the continuous problem described above is to find the 
unique positive flux Og), and the corresponding positive number k for 


which the system 1.1 has a solution subject to the conditions lee. In 
practice, we replace the continuous system by a finite -difference approxi- 
mation and solve the resulting system. To this end, let the reactor domain 
be defined by Rj = R= Rp, @_p = 9 = Oy. In addition, we lay down a mesh 
Ry, = Ro < R,<... <Ry = Rr; and Op = Oy 04 <.ee OS Om. Here we 
include the possibility of 9g = 0,9@q = 2m, which is the full periodic case. 
The mesh is assumed to be laid down in such a manner that region bound- 
aries occur only along mesh lines. Consider an interior cell centered at 
(r;, 85), as shown in Fig. 1. Integrate each member of the system 1.1 over 
this cell, apply Green's Theorem to the divergence term, and make the 
approximation 


f f g(r.) rdrd@ = $p(rj,6j) Vij = Pijg Vij: (1.3) 
3 ie 
where 
<a ; 2_ R24 
Vij = (A0j;/2)(R2-RE_,) 


is the volume of the cell and 
AQ; = 0; = os 


Having done this, we obtain the system 
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The surface integral extends over the four surfaces of the (i,j)th cell. 
Along the surfaces extending from points | to 2 and 3 to 4in Fig. 1, we 
have d¢/dn = d¢/dr. Region boundaries may lie along any of these sur- 
faces, so that ordinary interpolation will not suffice for even a first-order 
approximation. However, if we impose the conditions that the flux and 
current are continuous across the region boundaries and if we neglect 

the dependence of d¢/dr on @, we can approximate two of the surface 


integrals as follows:° 
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Here we have dropped the group index, and we are using the subscripts 0, 
E, and W for the subscripts (i,j), (it1,j) and (i-1,j), respectively. The 


linkage coefficients are defined by 


: ( L ) a eed geet 
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E e E 


(1.6) 


where Dy, Dp, and Dw denote the values of the diffusion lengths in the 


respective cells. 


Fig, 1 Interior Cell, Centered at (r1j,6)j) 


For the remaining two sur- 
faces, d¢/dn = (1/r)(d¢/90), and 
an approximation is sought for this 
tangential derivative along the two 
surfaces. Consider the surface 
extending from points 2 to 3, as 
shown in Fig. 2. To approximate 
the tangential derivative at the 
region boundary, we first approxi- 
mate the derivative of ¢ in the 
direction of the chord joining the 
points 0 and N. This is done just 
as before, where now 
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A 
This approximation (D/L) (dx - $0), is then projected in the direction of 
the tangential derivative and the resulting component, 


A A@:4, - Ad; 
(2) (On ~ 0) cos A, 
N 


is taken as the approximation to the tangential derivative. Thus, 


3 A 

BOes = Des 

f D = dQ = (2) (On - 9) AR; cos ee. (1.8) 
2 


If we set dq = Oo + dy then 


Ag; + AGi 4, 
ou = (R; + Runs) sin a7 ae ; (1.9) 
and 
b47 6 
UU 8 
In = Se ae RE (1.10) 
N 148, 1+, 
where 
sin AG, 4, 


= 
H 


U : ° 
sin 3 AG; 


The surface integral over the surface from 3 to 4 is handled in an analogous 
manner. Having finished the case of interior mesh cells, we must now con- 
sider the case when one or more of the surfaces of the cell are part of the 
exterior boundary when a boundary condition of the form A d¢/dn + Bd = C 
is imposed. If, for example, the surface from |] to 2is an exterior surface, 
we obtain in the usual fashion® 


2 

~ D(C -B¢ : 

po an ~pat ak F (1.11) 
1 A +i AR;B 


with an analogous expression for the integral from 3 to 4. If, for example, 
the surface from 2 to 3 is a part of the exterior surface, then AG ia, = 
therefore, Jn = 0, oy = 0%; and W®,y = 0. Thus, 


3 
D(C - B AR; 
poe dQ = Poe . Peek eéae AG: (1.12) 
; on ae, 


A + 69B 


where 
- a | 
So =~ (R; + R;_)) sin 4 AG.. 


Note that if an inhomogeneous boundary condition (C # 0) is present, then a 
source term is present and the relevant calculation is a source calculation. 


Fig. 2 
Detail of Interface between Two Interior Cells 


Having approximated the four surface integrals for the case of both 
interior and boundary points, we can make the following definitions for an 
interior cell: 


(2 
aij = a Ri Ae), 
NN 
(2 \ AG; > AG; 
bi; = \—] OR; cos ’ : 
4) NE IN 2 (1.13) 
oh meas ae 
and 
Si Paget 


At the boundary cells, we make the following definitions: 


DBR, A?. 
ara = ee Ia ea 
‘J A++ ARB 


DBAR, 
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DBR, AG; (1.14) 
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and 
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In the full periodic case, 6) = 0 and OF = 27; i.e., the top and 
bottom boundaries coincide. Thus, in this case we have 
d;, = bis lsisl. (1.15) 


Here we have dropped the group and region index on the quantities D, A, 
and B. Finally, if we set 


_ = as: o ag - 8 2\v.. 
Cijg = Fijg + Bile + Cijeg + dijp + (08 + DB ) Vij 
the finite-difference system 1.4 takes the form 


-(aijoditt, j,g a Dijgdi,jti,eg + Cijpi-1,5,2 di jg¢i,j-1,¢) 
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Here, since we have assumed that C = 0, the calculation is a reactivity 


calculation rather than a source calculation. The case of a source cal- 
culation will be discussed in Section IV. 


B. Finite-difference Equations XY Geometry 


For completeness, we shall exhibit the expressions for the finite- 
difference coefficients for XY geometry. Let the reactor domain be 
defined by 


Again we construct a mesh 
Xr = Xo < Xi < eS Xy = Xp; Yp = 2 ne oes a <Yy = Yr 
in such a way that region boundaries lie along mesh lines. Letting 
Ax, = Xx; - an and AY; = Yj - Yjey 
we set 
ais oP ijg Pitt, j; gs 
“ijg ~ 2X;D; 


et Oh Die 


b a ee (1.17) 
Ne AY Di jig + AY 4 Daj | 
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= Ai-1,j,2° and dijg = bi j-1¢2 


Cut 
1jg 


at interior mesh cells. At the boundaries, we set 


Dee bs. ‘ 

ae ee 7) 1<j 3: 

I Se Oe ee ee ee oe = —™ 9 
WB Any +2 AX Bye 

ee = ea eee 

®  Aisg +3 AYsBisg 

(1.18) 
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1Jg 


and 
_ DigBirgAXi 
ilg 


A + 5 AY,B; 


i1g 11g 
where the boundary constants Aijg and Bijg are permitted to vary by region 
and group. The form of the resulting finite-difference equation is the same 


as that in r@ or rz geometry. 
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Il. ADJOINT CALCULATION -._ 


The differential operator associated with the partial-differential 
system 1.1 and 1.2, when viewed as an operator in the space L* over the 
reactor, has an adjoint operator in the same space. If the boundary con- 
ditions 1.2 are homogeneous, then the domain of the adjoint operator 
coincides with that of the real operator. Thus the functions that lie in 
the domain of the adjoint operator satisfy the same boundary conditions 
as the functions that lie in the domain of the real operator. Hence, the 
adjoint operator is, in this case, just the formal adjoint operator. Con- 
sequently the finite-difference approximation to the adjoint can be found 
by forming the transpose of the real finite-difference equations. Before 
displaying the transposed system, we shall put the finite-difference equa- 
tions in the form of a matrix equation. The ordering of the space will be 
that of channel ordering,’ which is just a permutation of the usual group 
ordering. The ordering is defined in the following way: For eachi, 

1 <i <I, the set of points 1, = {(i,j,g): 1 = j<=J, 15g = Ghis called 

a channel. The points within a channel are ordered first with respect to 
the energy groups; then within each energy group they are ordered with 
respect to the points in the channel. Thus the flux vector @$ is partitioned 
in the following manner: 


fo) cae Dos cers oz), 
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Oe = Ais Disses OIG) l<isI; O51). 


Pig = (Dior Pipes soy Pigg) }S=i5 I; 15 g = Cr 


where the prime denotes the transpose, since all vectors in this discussion 
are column vectors. Relative to this ordering of the flux space, we shall 


define the following matrices: 
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(2.2) 


(In the full periodic case, Jig has the element -bigJ appearing in the extreme 
upper right- and lower left-hand corners.) 


Kig = diag[aigs. Aigo)? cues aigs|3 


t t t 
Bee - diag|V, Be aka vizofs" |) 


11°11 12°12 


br} 
7 
GQ 


& 
i a diag |V il oF, 8 (vo,)8 S09 V; 7085 *(vo,)§ a 


Then, collecting the above matrices, we define 


21 Gl 
Jin -B; ad -B; 
Bi’ Ji2 
a J5 = , 
1G 
-Bi Jig 


b. Kj = diag [Kis Kia: eeed Kic], 


and 
11 21 G1 
Fi Bi Fi 
12 22 
EF Fy 
res BF; = 
1G GG 
Be ee Re ener ree reer eee FE; 


Finally, we define the following matrices: 


Jy -~Ky 


(2.53) 


(2.4) 


(2.5) 
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In terms of the matrices M and F, the finite-difference system 1.16 can 
be put in the following form: 


Mo = - FO. | (2.6) 


As we have noted (see p. 14), the assumption of homogeneous boundary con- 
ditions enables us to approximate the continuous adjoint by simply solving 
the adjoint equations of Eq. 2.6. Thus we define o* by the following equation: 


M*o* = i FXO, (2.7) 


i 


where the * denotes the operation of taking the transpose. Forming the 
transpose of M involves taking the transpose of each of its blocks, and, in 
particular, would involve the transpose of the matrice J;. In the great 
majority of cases, J; is block lower triangular, and the iterative algorithms 
are designed to exploit this fact. However, J will then be block upper tri- 
angular, which is undesirable if we wish to use exactly the same algorithms 
in both cases. To circumvent this difficulty, we form the adjoint matrices 
and then invert or reverse the order of the energy groups. Hence we define 


the following matrices: 


as Jig = Ji G-gt.? 
b. Kig = Ky G-gta 
(2.8) 
— t ! 
C. Bee a Bes pre aet 
4; ee oe eee, 
Then with es F., K, M, and F defined in terms of these matrices as 
before, we can solve the equation 
= i 
My = = FY, (2.9) 
which has the same form as the real equation. 
We obtain the adjoint flux $¢* by setting 
(210) 


* - 
%i 08, = Vi G-g41,j° 


IIl SEARCH STRATEGY AND a-CALCULATIONS 


Consider a critical system; at time t = 0, let the system be modified 
in some manner. Suppose that this modification or change is a constant in 
time; then if the matrices M and F reflect the conditions of the modified 
system, the flux ®(x,t) can be described by the equation 


Vol 2 (ee) Oo, (3.1) 


where V is a diagonal matrix whose elements are the speeds associated 
with the energy groups. After a sufficient length of time, the flux can be 
reasonably described by a solution of Eq. 3.1 separable in time. Letting 
O(x,t) = (x)e% in 3.1 leads to 


(M+av7!)@ = F@. (3.2) 


Note that since aV~! is a diagonal matrix, its effect is to change the diagonal 


elements of M from e,,; to Cig] + Vy a Vij Thus the form of Eqs. 3.2 is 


identical to that of Eq. 5 6 eit " = ae if we call the result of solving 
2.6 a k-calculation, we see that a means for solving 3.2 can be based ona 
series of k-calculations. Thus, for each given a, we solve the system 


(M+av7!) o(a) = os Fo(a). (3.3) 


Then we search for that value of a for which k(a) = 1. Thus we see that 
an a-calculation is just a particular type of criticality search. 


To discuss the strategy employed in a search calculation, we shall 
have to describe briefly the method used in solving Eq. 2.6. Now recall that 


eae ee ee ee aioe) D.._B 
Sig)? “aie. weit igi t igs * Yu (oR ij ijg | 


where for each (i,j), we have 


t 
oh = of fs Sof 
g'#g 


Thus, 
Cig} > bigi + dijg = Pigj + Pigj-1: 


| : ; eS 
hence, J: exists and is positive. Moreover, 


1g 


1? 
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hence, sz? exists and is nonnegative. As a consequence, the matrix M7 F 
exists and is nonnegative. The problem posed by Eq. 2.6 can be viewed in 
the following manner: Find the largest positive eigenvalue and the cor- 
responding positive eigenvector of the matrix M7!F. The method used to 
find these quantities is that of power iterations.* The method consists of 
the following iterations: Starting with any nonnegative vector @ °) and any 
positive number k °), the sequences 7%"/ and k\"™/ are generated as follows: 


) (ya) yl) (3,4) 


n) — etl = (n-1 n) - n-1 
y\ = (M EF) y ), kc a Kc (y(n) y(n-1)) ) 


k(n) 


where 


(1.9) = D) Vigitigi- 
ijg 


Then, if the matrix M7'F is nonnegative and irreducible, it is well known" 
that the sequence k\"/ will tend to the unique largest positive eigenvalue of 
the matrix, and the sequence 7~/ will tend to the corresponding positive 
eigenvector. In the above situation, we see that the generation of y\”’ from 
y(n-1) involves the solution of the equation 


1 
k n-] 


Myln) - pyln-1), (3.5) 


In general, the matrix M is far too large to be held in fast core; however, 
it is block tridiagonal with nonpositive off-diagonal blocks and nonsingular 
diagonal blocks, which are diagonally dominant and have nonnegative in- 
verses. Thus, a Gauss-Seidel iterative procedure will converge toa 
solution of Eq. 3.5. Referring to Eq. 2.5a, we define 


a. D 


il 


diag|J,, Jo, ooed J4Js 


O O 
Ky OR b (3.6) 


b. Ge © - U = 
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Then M = D-L- U, anda Gauss-Seidel iteration relative to this splitting 
is defined by the following procedure: 


(D- 1) x4) = ux-1) FuUM-1) gf = 1,2, ..., (3.7) 


is 1 
(2-1) 


with x0) = yln-1), The sequence x4) will then converge to the solution 
w\n) of Eq. 3.5. Moreover, if there is no up-scattering, the matrices J; 
can be inverted directly since they are then block lower triangular, with 
tridiagonal matrices forming the diagonal blocks. If up-scattering is 

present, the matrices J; satisfy the necessary criteria so that a Gauss- 


Seidel iterative procedure can be applied to these matrices. 


The present code does not go beyond the first iterate in the iterative 
procedure defined by Eq. 3.7. Thus, y(nt 1) is defined to be equal to x0), 
hence, the iterative procedure used in this code is as follows: Starting with 
any positive number k °/} and any nonnegative vector °), a sequence of 
numbers k\”/ and vectors ¢\"/ is generated by 


) (g(n) g(x) 


(g(n) g(n-1)) 
(3.8) 


(D = 15) gin) = ugln-1) ia ! ro(n-1) | (0) = j(n-1 
| (n-1) 


Thus we see that the actual iterative procedure used in this code is a non- 
stationary, nonhomogeneous power iteration, where the product 


n 
TT Sy 
L=1 


replaces 


Here Sy denotes the nonnegative matrix (De iy" (0 +k) F!: 


Before we discuss the effect of the above iterative procedure on the 
search strategy, we shall display the above iterations in more detail. Thus, 
using the definitions of D, L, and U, we have the sweepthroughthe channels. 


A) a K.g(27}) = F.g(h-}) fori =1,2,...,1 (3.9) 
k 
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If up-scattering is present, the matrices J; will not be block lower 
triangular; thus, the vector @s$"/ will be found iteratively. Within a given 
channel, we have the following iterative scheme: 


t t 
evr) = DY pe eyim), Do BP eyim)s Ky, itl 


gs g-g' | 
(n-1) , 1 : g'g.(n-1) 
n-1 n-1 = 
so Ki pPitig + (n) > FE Pig for g = L; Ze ee G. 
k gi=l 
(3.10) 

Here we take yo) = gin-t) and upon achieving convergence, we set 
g(a) = y(m) for ¢ = 1, 2,255 G, and then move on to the next channel. 


Pho up- scattering is present, we obtain Oe exactly in one sweep through 
the groups. Except for the full periodic case, the vectors YM? are found 
by a direct inversion procedure based on the factorization ofa tridiagonal 
matrix into the product of a lower bidiagonal and an upper bidiagonal 
matrix.!’? This technique is sometimes known as Choleski's method.’ In 
the full periodic case, the solution is found by means of an iterative tech- 
nique, which will be discussed in Section VI. 


Having thus displayed the iteration procedure at various stages of 
detail, we shall consider the question of how a search is carried out. In 
a search, we seek the value of a parameter such that the resulting system 
has a kere equal to a prescribed value. Thus, we seek that value X) such 
that k(X)) = ky where ky is given. This is accomplished by means of linear 


interpolation; hence, given 
(x(5),k(x(8))) and (xx(8-1),1(x(s-9)), 


a straight line is passed through these points, and x(s+1) is taken to be the 
abscissa of the intersection of this straight line with the line k = kp. The 
major problem is that an estimate for k(X) is expensive in machine time; 
moreover, the cost is higher than usual in this code for the following 
reason. Experience has shown that generally the rate of convergence of 

a power iteration procedure based on Eq. 3.4 is relatively fast (r = 0.9) 
for the flux and even faster for the sequence k\"/, (In this context, by a 
rate of convergence r, we mean that the error eventually behaves like 

r™, where n is the iteration count.) Since the matrix M is so large, we 
cannot perform an iterative procedure based on Eq. 3.4. On the other 
hand, we could base a practical iterative procedure on Eq. 3.5 in conjunc- 
tion with Eq. 3.6. Then the same rate of convergence would apply tokq. 3.95 
however, the rate of convergence of Eq. 3.7 would be very much slower. 
Since the determination of the sequence k\"/ is based on Eq. 3.5, we would 
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expect the rate of convergence of (n) to be comparable to a procedure based 
on Eq. 3.4, provided a sufficient number of iterations are performed in 

Eq. 3.7. However, since the flux is determined by Eq. 3.7, we would expect 
its rate of convergence to be slow. Thus, insofar as rates are concerned, 
we would expect the sequence k n) to converge more rapidly than the fluxes. 
The present code uses only one iteration of Eq. 3.7 per iteration Of Eq. 3.23; 
thereby creating the iterative procedure, Eq. 3.8. In this procedure, the 
rate for the sequence k\"/ is now comparable to that for the flux estimates, 
hence converging at a rate of 0.99 or greater for the average-sized problem 
(Appendix A, Sample Problem 1). This situation, namely the slow conver- 
gence of the estimates k\"/, causes an extreme amount of difficulty in the 
search procedure by compounding or accentuating the following dilemma: 
On the one hand, for a given X 1), a sufficiently accurate estimate k(n) (x (i)) 
must be found for «(x (4) to ensure that the sequence x (i) will converge 
to Xj. On the other hand, because of the curvature of k(X), an accurate 
determination of k(x 1) does not lead to a comparably accurate estimate 
x(it1) of Xp when X 1) ig not near Xp. Taking these considerations into 
account, we have developed a search procedure based on the following 
observations: 


1. Each (2) (3) is very costly; its value in time corresponds to a 
complete sweep of the mesh for all groups. 


2. The sequence k(n) (x) is slowly converging, its rate being of the 
same order as the flux. 


3. The sequence (2) (x) changes slowly as a function of X. Thus, 
if «(n)(x) is an acceptable estimate for a given X and if X' is the next 
estimate for X,, then k\'/(X') is close to k\”/(X). 


4, After a sufficient number of iterations, the sequence «(n)(x) 
converges to k(X) in a geometric fashion. 


With these observations in mind, the following strategy is employed: 
Suppose that we are beginning the iterations needed to approximate ki = 


(x (i)) corresponding to a control value x(t), Here x (i) may be an initial 
guess or the resylt of interpolation. The code then generates successive 
estimates ky : ki? , ..., kS]) for k.. Because of observations 2, 3, and 4, 
no change is contemplated until at least 15 iterations have been completed 
and a minimal level of convergence has been achieved. We then determine 
whether the sequence ki is heading towards the desired value k). If the 
sequence is heading away from kp, we stop the iteration as soon as the 
following criterion is satisfied: 


icf) - ee] < 0.2 [abr) - & I. (3.11) 


Zz 


Of course we do not know kj; therefore in) ~ k, | must be estimated. Based 
on the assumption that the convergence of kim is geometric, this difference 


is approximated by 


(n) Ss l n) 
eer 2 ae py aae (3.12) 
1 

where 

»(n) ae lay 

i Jawa) 

i 

and 


The test defined by 3.11 in conjunction with the approximation 3.12 is made 
only after at least three successive unextrapolated iterations have occurred. 
This restriction is necessary because extrapolation destroys the geometric 
character of the convergence. The case in which the sequence k\"/ is heading 
towards ky is somewhat more complicated. In this case, an effort is made 
as soon as possible to determine whether the iterates k\?7 will pass through 
the desired value ky. Again this determination is made under the assumption 
that the sequence in) 3 is converging geometrically; thus, the same restric- 
tion with regard to ‘flux extrapolation is in effect. If the se ks 1} for k, 
lies on the other side of ky or is close to kg (that is, Ik, (n) < (10)e os 
then the code continues to iterate without making a Sontnel ‘change. On the 
other hand, if the estimate k}"/ lies on the same side of kp as ki? and is not 
close to ky, then the test defined by 3.11 and 3.12 is applied. If the test is 
satisfied, a control change is made using the estimate ki") for k.. If the 

test is not satisfied, then the code continues to iterate (see Appendix A, 


Sample Problem 3). 


Having discussed the decision-making process involved in seeking 
a value of the control parameter (X) such that k(X) = kyo, we turn next to 
the logical control of the search process. Of course, what we ultimately 
seek as the result of a criticality search is not the required order of X but. 
the value of some reactor parameter(s) such that k (parameter(s)) = kg; 
i.e., we must have parameter(s) = F(X). The functional relationship used 


in the program is 


pis) - p(y ‘ x(S)sp_), 
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where 

X° is the value of the control parameter for the sth pass, 

pl?) are the initial values of the reactor parameter(s), 

p(s) are the values of the reactor parameter(s) for the sth pass, 
and 


5P,, are the parameter(s) modifier(s). 


The code currently allows the following four choices of the parameter(s) 
P,, (i.e., four different types of criticality searches are possible): 


1. Composition Search 
Py = VEm,c: 
where VFym,c are the volume fractions for materials m in composition c. 
2.. Dimension Search 
Py, = AX; and/or AY;, 


where AX,(AY);) are the mesh increments (i.e., interval lengths) in the 
X(Y) direction. 


3. Buckling (B*) Search 
Pn = BY, 
where B¢ is the transverse buckling for region r. 
4, a Calculation 
Pay SO; 
where a is the asymptotic inverse reactor time period. 


The logical control of the criticality search is essentially as given 
in the CANDID2D document,! but is duplicated here for completeness. 


A first guess, Xa; and either a second guess, Xp, or an estimate 
of dk/dX, (KDOT), and ky are required input. If KDOT is given, then Xp 
is computed from the following equation: 
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Xp= Xat (x, - k(x) /KDOT. 


Once k(X,) and k(Xp) are computed, linear extrapolation is done until 
ky is bracketed. A maximum number of extrapolations to bracket Ky is 
provided in the input. Once ky has been bracketed, the values of k(X) 

and X that best bracket ky are stored in Xy» ky» XR; and kp such that 


be hy 1 We 


For the purpose of the following discussion on bracketing k, let 
us assume that X < Xa < XB < X. If an extrapolation results in X < X, 
then a midpoint is attempted for X from A = (X + Xa)/2. If a later ex- 
trapolation has the same result, X < Xx, then the bound itself is used; 
126 e Kw. SK, : 


A third failure at the same bound will be a terminal error. The 
above discussion follows through in a similar manner if XB < XA and/or 


the failure is at the other bound, X. 


X,, kj, X2, and k, are the interpolation parameters and are chosen 
in the following manner: Once ko has been bracketed, set A, and X, equal 
to Xj, and XR, and interpolate for X from the following equation: 


Ko xs k(X)) 


K(X) te x) (X, - X)); 


X = X, + 


then k(X) is computed, and the bounds are updated as follows: 


If k(X) < kp, then Xj, = X and ky, = k(X); 


or 
if k(X) > kp, then Xp = X and kp = k(X). 
The closest two values of X are used for the next interpolation: 
if |x - x,| < |x - x,|, then X, = X and k, = k(X); 

or 


if [io RG |S | then XS X-and: kp. =: K(X). 


Now interpolate for a new value of X. If X does not fall in range 
of the best known bounds, 


<< Xpand XR 
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or 
xX > Xz,, and Xp, 
then the bounds are used: 
X, = Xu, k, = KL, 
or 
X, = Xp, k2 = KR; 
this will yield a valid X. 


This procedure is continued until either the convergence criterion 
is met or a specifiable number of interpolations is exceeded. 


2.6 


IV. SOURCE CALCULATION 


A source calculation can arise in two ways. Either a fixed source 
is prescribed, or an inhomogeneous boundary condition is present. In 
either event, the basic assumption of a source calculation is that the reac- 
tor system in the absence of the fixed source is subcritical. In terms of 
the matrices introduced in Section II, a source problem can be written in 
the form 


(MeaF).o-= 4, (4.1) 


where f is a prescribed nonnegative source. The assumption of subcriti- 
cality ensures that the matrix M - F is not singular. In addition, it ensures 
that the solution ¢ is nonnegative when f is nonnegative. This last state- 
ment is by virtue of the fact that 9 = Clie MIF)? M'£, and by the hypothesis 


that the spectral radius of M-'F is less than one. Thus (1-M7’F) 2 0. 


The method used to solve Eq. 4.1 takes advantage of the existing 
method for solving the homogeneous problem. Hence, in terms of the ma- 
trices introduced in Section III, the iterative procedure employed in the 
absence of up-scattering can be written in the form 


(D-L) om) = (U+F) gin-1) + f. (4.2) 


As far as the mechanics are concerned, this procedure differs from 
the kerr case only to the extent that the eigenvalue estimates k(") are not 
used in the fission source. Recalling the structure of D - L from Section Ill 
for no up-scattering, we see that (D- 1)! exists and is nonnegative. Fur- 
thermore, U + F is a nonnegative matrix; thus, the splitting M - F = [DetLa< 
(U+F) is a regular splitting.” As a consequence, the spectral radius of 
(D- L)-(U+F) is less than one; thus, 4.2 is a convergent iterative procedure. 


When up-scattering is present, we shall assume that a fixed number 
of up-scattering iterations are performed. From Eq. 3.10, we recall that 
within each channel i, 1 =i < I, the up-scattering iterations have the 


following form: 


(m 'g._(m) Joe cooles 6 ames (n) 
gi<g g'>g 


G 
(n-1 lige (m1) 
z3 Kigdt tig + y FP Ebi g! +. fig g = 1s Z, ares G. 
en (4.3) 


Here, 


(o) _ 


(n- 1) 
1g , 


and at the end of £ iterations we set 


Referring to the definition of Jj in Eq. 2.4a, we form the block decomposi- 
tion Jj = Dj - Ly - Uj, where 


D; = diaglJi,, 555.0889 Jicl, 
and L; and Uj are the remaining block lower and upper triangular matrices, 


respectively. In terms of this decomposition, the up-scattering iterations 
take the form 


(D; - Ly) vf") = uy) 4 Kyo '8) + Kyo) + Fg) 4 (4.4) 
If we set R; = (Dj - L;)7' Uj, and recall that y°) = oo) and of) = yl"), 
then oi?) has the following form: 

ae peg) 4 (r- Rf) 5724 Ki 10 + KO) 4 Fgh) 4 af. (4.5) 


We can write this iterative procedure in a manner analogous to Eq. 4.2 if 
we define the following matrices. Let M and N be defined by 


I 
i SRe) as I é 


M = (4.6) 


£ 
C) roe L =I ae 
-(I- Ry) Jp Ki-1 I 


and 


Rae) 3 (1 RY) 37K, C) 


w 
N = rR? + (1-R#) 53°F, (I- R4) 32"K - (4,7) 


~~ 


& RO -87) 0 Py 
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With these matrices defined, the iterative procedure 4.5 can be written in 


the form 

mo) = No) + g, (4.8) 
where 

gi = (1-R?) 3738 oy oie esa ae emer 


At this point we can make several observations: 


1. For each channel i, the matrix (1 - Rr?) te is nonnegative and 
nonsingular. 
2. The iterative procedure 4.8 is convergent. 


3. The solution thus obtained satisfies Eq. 4.1. 


Starting with the first observation, recall that Rj = (Deeg) is 
the block Gauss-Seidel iteration matrix associated with the matrix Jj. The 
matrix J; has on its block diagonal nonsingular matrices which have posi- 
tive inverses. Moreover, the block off-diagonal matrices are nonpositive; 
thus, R; is nonnegative and has a spectral radius less than one. For the 
Same reasons, Jz? and (Dj- L;)"' are nonnegative. As a consequence, we 
see that Jj! = R:Jj' since 


(I- Ri) JZ? = (I- Ry)(I- Ry)" (Dy - Ly)! = (Dy- Ly)" = 0, (4.9) 

From this it follows that (1- rR?) J;' 20 for any £21. Since the 
spectral radius of R; is less than one, the same is true of R:; thus,{I- Ro) 
exists and is nonnegative. Hence, \I- Rj J; is not singular. 


Addressing ourselves to the second observation, we see from the 
first observation that M7‘ exists and is nonnegative; moreover N is also 
nonnegative. Thus the matrices M and N used in 4.8 form a regular 
splitting’ of the matrix M - N; consequently, the process is convergent. 

Concerning the third observation, suppose that yw satisfies 

A A 
(M-N) y = g. (4.10) 


AN N 
Now M - N = SJ7}(M-N), where 


J S25 diae| Tides sgl 


and (4.11) 


From the first observation, the matrix SJ~! is not singular; moreover, 
g = SJ-4f. Hence, w satisfies Eq. 4.1; therefore w= ¢, where @ satisfies 
the equation 


(M-F) @ = f. (4.12) 


Based on the above observations, we can draw the following conclu- 
sion: When up-scattering is present, the solution does not depend on the 
number of up-scattering iterations performed. The use of up-scattering 
iterations would then rest on the possibility that they will improve the rate 
of convergence or at least the asymptotic rate. The possibility of improve- 
ment rests heavily on the behavior of N defined in 4.7 as a function of the 
number of up-scattering iterations, 4. However, although the asymptotic 
rate of the process defined by Eq. 4.8 does increase with increasing 4, it 
is still not sufficient to offset the consequent increase in computation time. 
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V. ACCELERATION 


In this section we shall describe the acceleration scheme that was 
implemented for source calculations. The scheme used was two-term 
Chebyshev extrapolation applied to the iterative scheme 4.2 (or 4.8 for 
up-scattering). This scheme was chosen because a two-term Chebyshev 
procedure is used to accelerate the kes calculation. Thus a minimum 
amount of effort was involved in adapting this procedure to a source 


calculation. 


Recall from Section IV that the iterative procedure used to solve a 
source problem can be written in the form 


g(nt1) - pgm) 4p, (5.1) 


Here R = (D- L)"(U+F) and b = (D- L)"! f, if no up-scattering is present. 

If up-scattering is present and we assume that a fixed number of up- 
scattering iterations are performed, then R = MIN and b = M ‘f, where 

M and N are defined by Eqs. 4.6 and 4.7. Neither the matrix R nor its 
associated Jacobi matrix is symmetric; thus, to justify the use of Chebyshev 
extrapolation, we shall have to assume that the matrix R has real eigen- 
values. The two-term Chebyshev extrapolation* can be described briefly as 
follows: Let @ be the solution of Eq. 5.1, and let e(™) = ¢ - g("), It then 
follows that 


e(p) = Re(P-}) - RPe(°), (5.2) 


If we assume for simplicity that the eigenvalues X4g of R are simple with 
1 >A); >A2>... 2 O, then e(P can be expressed in the form 


ai) = e,rPu, t corbu, +... (5.3) 


If, on the other hand, we had used 


ied R-p.) R - 
yr 1 - Ap 1 - Ap.) J 04 


in place of RP on e (0), we would have obtained 


j od 
~~ + C2 l] a ong vee 5 4ud.cs (5.4) 


Here we can See the possibility of making e(p) as defined by Eq. 5.4 smaller 
than e(P) as defined by Eq. 5.3. For if we set 


Pee a0 
Qp_(A 5 Aa) a TT T. —, (5.5) 
j= 


and then choose Q_A; X,) to be the minimal polynomial in the supremum 

norm of degree p over the interval [0, A,], then e(p) as defined by Eq. 5.4 
will be smaller than e(P) as defined by Eq. 5.3. The solution to the above 
minimization problem is given by 


2r 
op(2- 1) 
oF) (5.6) 


where C,(x) is the Chebyshev polynomial of degree p normalized to be ! at 
x = 1,. The roots of the polynomial in Eq. 5.6 are then given by 


oP) = ab + cos {20-9 + 1) zt for 4 SUZ, 4445:'p; (5.7) 


Since the matrix polynomial Qy(R; 1) is expressed in its factorial form, we 
accomplish the extrapolation by applying one factor at a time. Suppose that 
starting with a flux ¢ n) and an estimate De for A,, we decide to perform an 
extrapolation cycle of order p. The roots {ofP)}P are then computed from 


Eq. 5.7 by using oe in place of X,. For use in the actual application, we de- 
‘Then we generate in succession 


fine the numbers pip) = ( - o y : 
Ai) 2 aeeads 


fini) - —1 (glnts) - a(elgla)) = plrlelat) + (1 - ofp) ol), 


gin+2) _ Rint?) +b: 


*(n+2 | n+2 (n+ p n+2 p A(n+] ( : ) 
aN ) = ( 16 Ve ale) ol )\ ph pl 4 ( 2p\ ) an ), 

ran ] At 3 n A n a4 

g(ntp) _ eRe (s(n) zs o(P)p(atp ) as piP)g\ +p) 4: (1 é a(P)) al ee ) 


A 
With the generation of ¢(n+p) we have completed an extrapolation cycle of 
order p. The error vector, @(ntp) - d - é(ntp), then satisfies 


A(ntp) - Qp(Ri X1) ain) (5.9) 
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As is well known,*’* the effectiveness of this extrapolation is sensi- 
tive to the accuracy of the estimate ,; for X,;. In this code, the initial esti- 
mate for ,; is found from the fact that, for sufficiently large n, 


jg(nt1) - gl) |], | 


i 5.10 
}o(2) - gla-a)|], wan 


]> 


where 


Isle = ( 2 Oh) 


ijg 


Subsequent estimates are made using an up-dating procedure introduced by 
Varga.” This up-dating procedure is used in the following way: Suppose an 
extrapolation cycle of order p was begun with the nth iteration, so that 
g(ntp) is the result of this cycle. An unextrapolated pass is then made gen- 
erating g(ntpt!), and the following quantity is formed: 


|| o(ntp+1) - ¢(ntp)| e 


E = —ooeeeeeeeeneneneeo 5.11 
mp = g@H)_ gp, ie 
It can be shown? that 
En,p = |Q_(% 15 Aa) | (5.12) 
Now 
a, —ae | 
Q,0; ri)| = Q,(0i5 41) = azo ySCtiCéiforr OA = A, (5.13) 
Cp (= - ! 
ry 
and 
Q,(A5 Aa) = Q,(Ai; Ai) for ve oar (5.14) 
Hence, using Eq. 5.12, if 
1 
E (5.15) 


—— 
Se 
PX, 


then we can conclude that A, = oe and also that the extrapolation was as 
effective as we can expect since 1/[C (2/r, -1)] is the maximum of the 
modulus of QA; 1) in the interval lO, ee: On the other hand, if 


ES 


] 
C : | 
PLA, 
then, using Eq. 5.14, we can conclude that ne < A,. In this case, the esti- = 
mate for i, will need to be moved up. Since dX, < dj, we know that Q”(Aj; r,) 


is positive; thus, using Eq. 5.12, we define the new estimate KG as that value 
of X% for which 


Baa (5.16) 


Page O5(ls ay): (5.17) 


This yields a new estimate a for which a < pe <i; 


A source calculation then proceeds in the following manner: Starting 
with an initial guess for the flux, we perform 15 unextrapolated iterations, at 
which point we make an initial estimate 1 for A, using Eq. 5.10. If this 
estimate, which is also a measure of the rate of convergence, is greater than 
0.75 and less than 0.995, then we perform an extrapolation cycle of order 6. 
If the estimate is greater than 0.995 but not greater than 0.999, the order is 
limited to 5. On the other hand, if the estimate exceeds 0.999, then this 
value is used in place of the estimate. On the other hand, if the estimate is 
less than 0.75, then the unextrapolated iterations are continued until the 
estimates exceed 0.75 or convergence is attained. At the end of an extrapo- 
lation cycle, we make a test according to inequality 5.15 as to whether the 
cycle was effective. If it was not, we make a new estimate according to 
Fiq. 5.17. Three unextrapolated iterations are then performed before enter- 
ing the next extrapolated cycle. Generally speaking then, the pattern is 
extrapolated cycles of length p, interspersed with three unextrapolated 
iterations. Of course, there are exceptions to this general pattern. For 
example, if the rate of convergence (as estimated by Eq. 5.10) drops below 
0.75 or rises above one during unextrapolated iterations, then these unex- 
trapolated iterations are continued. 


The problem is continued until the following criteria are satisfied: 


[eI] "IIo - gina) < €% (5.18) 
and 
k(n) - k(n) < ey, (5.19) 
where 
(n) (n) 
git gi? 
™ ijg. a tack. de 
k(n) = max Gea) k(n) <= eS ae) (5.20) 
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and the extrema are taken over the set of points (i, j, g) for which ane) = Eg 
(see Appendix A, Sample Problem 4), Also note the following fact: a 
If any component of the flux should become negative during, or at the con- 
clusion of, an extrapolation cycle, this component is set equal to zero. The 
intuitive justification of this is as follows: The external source is nonnega- 
tive: thus it might be expected to be strongly biased in the direction of the 
dominant eigenvector of R, which is positive. Thus the desired solution is 
positive and thereby biased in the direction of the dominant eigenvector. 
Thus, if a component becomes negative in the course of extrapolation, it 
indicates a buildup of a lower harmonic. This means the approximation is 
heading away from the solution and will have to pass through zero on its 

way toward the solution. Thus, setting the negative component equal to 

zero can do no harm and indeed prevents the difficulty that arises when 


any component of the flux goes negative. 


Many of the procedures used to accelerate a source calculation have 
been incorporated into the Keff acceleration procedures. To describe these 
procedures, as wellas others that were developed, we shall describe the 
acceleration procedure in some detail. For simplicity in description, we 
shall assume that no up-scattering is present and that the problem is not a 


full periodic problem. 


Recall from Section III that a Keff calculation involves finding the 
largest positive eigenvalue and corresponding eigenvector of a nonnegative 
matrix. We saw from Eq. 3.8 that this code solves the problem by relying 
entirely on power iterations (disregarding up-scattering and periodic 
iterations). Thus, from Eq. 3.8, the basic iterative procedure is given by 


(a-1) Ao), o{)) , (5.21) 


(n) — py ! (n- (n 7 
g(2) = (D-L) urea |¢ i) xin) = k (4a), 23) 


If we set S(k) = (D- by) os (1/k) F] with Sn.) = s(ic(a-2)) we can write 
Eq. 5.21 in the form 


g() ss sin!) = TT so, .(n) Ze 1.(n-1) (5.22) 


m=1 


It is not hard to see that if k is the dominant eigenvalue and ¢ the corres- 


ponding eigenvector which satisfy 


FQ, (5.23) 


then 


S(k) @ = . (5.24) 


Moreover, 1] is the dominant eigenvalue of the nonnegative matrix S(k). Thus 
if we know that the sequence k n) as generated in Eq. 5.22 converges to k, it 
follows that the sequence g(n) will converge to ¢. However, the coupling be- 
tween k(") and g{n) is so complex that it is difficult to show that the proce- 
dure defined by Eq. 5.22 will converge. 


Let us assume that the process defined by Eq. 5.22 is convergent. 
The application of Chebyshev acceleration to power iteration is predicated 
on two basic assumptions: 


1. The eigenvalues of the matrix are real. 


2. The eigenvalue estimates are converging at a faster rate than 
the eigenvector iterates. Moreover, the extrapolation is not 
applied until the eigenvalue iterates have almost converged. 


Suppose both assumptions are satisfied; so that k(n) = k. Then the 
process becomes an iterative process with a fixed matrix, S = S(k). Thus, 


gl(Zti) - sp(Z) for 2 = 2, (5.25) 


where S has its dominant eigenvalue equal to one. Let 6 be an estimate 
for the second largest eigenvalue 0 of S, whose eigenvalues lie in the 
interval [0,1]. Then a Chebyshev extrapolation cycle of order p applied 
to the process Eq. 5.25, beginning with 4 = n, would yield 


plotp) = Q,(s; 5) o), (5.26) 


a) 
Q (55) = oe) 


The function C (x) is, as before, the Chebyshev polynomial of degree 
p normalized to bec l atx =.1. If we expand g(n) in terms of the eigenvec- 
tors of S, we see that the g(ntp) has the form 


where 


(5.27) 


g(ntp) = a,Q,(1; 5) wy + a2Qp(%; 5) ut ..., (5.28) 


where g(n) = ¢,}U,; + coUz +... and 1 > 0, 2 0,2... = 0 are the eigenvalues of 
S. Now Q,(1; 6) =1 and Qp(x; 6) is minimal over (0, Oo]. Therefore we have 
pertoemicd an effective extrapolation depending on where our estimate O is 
relative to 0). 
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We have described the standard technique for applying Chebyshev 
extrapolation to eigenvalue problems. Accounts may be found in Refs. 3 
and 6. Note that here the dominant eigenvector estimates play the same 
role that the error vectors did in the source calculation. Based on the 
above analysis, we could use the up-dating method for finding new esti- 


mates for 0). 


When the two-term Chebyshev extrapolation is applied to a keff prob- 
lem, the results are disappointing (see Appendix A, Sample Problems I! and 
2). This is in contrast to the source-calculation case (see Appendix A, Sam- 
ple Problems 4 and 5). There the results are quite good; for example, the 
extrapolation cycles make their theoretical error-reduction criterion 
(En p = Cp(2/r- iy) more often than not. On the other hand, for a keff cal- 
culation, an extrapolation cycle will almost never make its theoretical cri- 
terion, frequently does little better than the same number of power iterations 
would have done, and occasionally even causes apparent divergence. Of 
course, we can always recover from this latter contingency by performing 


power iterations. 


We shall explain why extrapolations in the keff case behave as they 
do. We shall also describe what remedial steps were taken in the use of 
extrapolation to minimize the effects described above. 


The problem in the keff case is that assumption 2 is not satisfied. 
That assumption 1 is essentially satisfied is shown by the following em- 
pirical evidence: If k = 1, then the matrix S(k) is essentially the same as 
the matrix R which appeared in the source calculation. In fact, S(k) = 
ARs (ice 1)(D - L)"? F. The source extrapolations were very effective, 
indicating that R had essentially real eigenvalues. Thus we would expect 
that at least for k near one, S(k) would also have essentially real eigen- 
values. This argument is somewhat weak since S(k) can be viewed as a 


perturbation of R by a nonsymmetric matrix. 


On the other hand, the evidence that assumption 2 is not satisfied 
is overwhelming. In almost every case, the rate of convergence of the 
eigenvalue iterates k(n) is of the same order as the eigenvalue iterates. 
Thus, it is not feasible to wait until the sequence k\"/ has converged before 
extrapolating. Let us consider the effect of performing an extrapolation 
cycle before the eigenvalue iterates have converged. Two effects are pres- 
ent. First, the estimate k\"/ for k is not accurate, and second, the iterates 
k(n) are changing during the extrapolation cycle. Consider the first effect. 
Assume that we have an estimate wu = k(2) and an estimate 6 for the domi- 
nance ratio of S(u) = s(k()), i.e., if A,(u) > AA(u) = ... are the eigenvalues 
of S(u), then we have an estimate for A2(u)/Ay(u). Note that since L Zk, 
ri(u) # L =o Agtk). et g(n) be considered in terms of the eigenvectors of 


S/(u); i-e., 


g(n) = cju, + Czu,t ..., ay aa (i) (5.29) 


Then, as a result of an extrapolation cycle, we obtain 
A - a 
g(ntp) - €,Q,(Ai(u); G) wy + c2Q_(Aa(H); B) uz +... (5.30) 


Figure 3 is a graph of Qp (x; 6). 


Fig, 3. Graph of Qp(x; 0) 


Now i(p) Z# ], and since pu is in general not even close to k, we 
see that ,(y) need not be less than one; or if it is less than one, it need 
not be greater than 6. Thus we see that the "reduction factor" 


Qp(r2(u); 6) 


Q,0..(u); 


Q 


(5.31) 


a 


need not be less than the theoretical reduction factor [C,(2/c - |) )"*-re2 
gardless of how well the dominance ratio 6 is estimated. Indeed we can 
see that for A,(u) < G5, the reduction factor can be greater than one. Thus, 
we can begin to see why the extrapolations are so ineffective, and how the 
higher harmonics can be built up during a cycle when eGo se anges 


The problem of a changing estimate i (n) during the extrapolation 
cycle adds an additional degree of complexity by preventing us from 
applying a Chebyshev polynomial. 


Having examined the causes for the ineffectiveness of the extrapo- 
lation procedure, we find that two courses of action are open. First, change 
the basic iterative structure. Second, modify the technique of applying the 
extrapolation procedure. The first course of action requires scrapping the 
computational portion of the code and starting from the beginning. This is 
being done. As far as this code is concerned, the procedures for applying 
the extrapolation have been modified to at least prevent disaster. 


The following procedure is available, on option, to help prevent the 
occurrence of A,(u) < 6. The code has a built-in sequence of numbers 
0.8, 0.85, 0.9, 0.925, 0.95, 0.96, ..., 0.99, 0.991, ..., 0.995. The user 
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specifies one of these numbers, say Og. Then in the first extrapolation 
cycle, the code will not use any estimate for 6 that exceeds Oy. If 5; is the 
next larger number in the above sequence, then on the second extrapolation 
cycle the code will not use any estimate © which exceeds 0,. This process 
is continued until 0.995 is reached. The above procedure then has the effect 
of restricting the growth of 0 as the iterations proceed. 


In this connection, a second technique is available for avoiding the 
situation A,(u) > G. This involves running a coarse-mesh problem and then 
refining the mesh and using the results from the coarse problem. 


The up-dating method for finding estimates for the dominance ratio 0 
is used only when the sequence k n) has essentially converged. As we have 
noted before, this seldom occurs before the flux has converged. The reason 
the up-dating method is not used before this point is that it forces the esti- 
mates for o to be large very early in the iteration process. Most of the 
estimates for 0 are made by means of 


jo?) - gln-1)]] , 


On, = Jom) - g-2)]|, = 0 (5.32) 


This monitoring is performed during the three- to four-power iterations 
that follow each Chebyshev extrapolation cycle. 
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VI. SOLUTION OF THE FULL PERIODIC PROBLEM 


The full periodic case arises in r@ geometry when the reactor do- 
main is annular. Thus the flux must satisfy dglr, 0) = bg (r, 27), This 
condition leads to the following form for the matrices Jig: 


Sigi “Bigg “Pigs 
“Dion “iga ~bis2 
Jig = : 6.1) 
“Big J-1 
“bi og “big ,J-1 Cig,J 


We recall from Eq. 3.10 that for each channel i and group g we 
seek the solution of the equation 


eee (6.2) 


ig dg = ig 


1g 
where fig denotes the entire right-hand side of Eq. 3.10. In the nonfull 
periodic case, bj,yz = 0, and the Choleski algorithm is used, Because of 
memory limitations and programming difficulties, Eq. 6.2 is solved 
iteratively. 


To effectively describe the iterative scheme used to solve Eq, 6.2 
in the full periodic case, we shall briefly describe the Choleski algorithm 
as used in this code. In the remainder of this section, we shall drop the 
subscripts i and g. The basic idea of the method is to factor the matrix 
into the product of a lower triangular matrix and an upper triangular 
matrix (Gauss elimination). Thus, starting with a system Ax = g, where A 
is a tridiagonal matrix, we obtain the pair of equations 


10s = 2; 
} 6.3) 
Ux = ®. 


Here L is lower triangular, and U is upper triangular and has the form 


-, OC) 
l aye 


c= l “Ys (6.4) 
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A forward sweep yields the sequences {yj} and {@;}; the backward sweep 
yields the solution {x; }. Thus if 


C1 -b; 
-by €2 -b2 
A = ; (6.5) 
-b 
pri “p 
then 
by 
.; eae ee! es gi/er; 
€1 
b: 
J . 
fe es OE.) Suds Opaae Roo ks (6.6) 
J ej = bj Vj-1 
Sor b, O. | 
OD; - i," for j = as 34 woey P. 
* jr Yj-1 


The backward sweep that generates the solution is then given by 


x = @~; | 
p 19 (6.7) 
x, = ee Pot for j Sop = by p= 2ion5 


Note that the sequence {w:} depends on the vector g, whereas the sequence 
{1.} depends only on the matrix elements; thus the sequence Ly;) is calcu- 
lated only once and then stored. 


We shall now consider the iterative method used to solve the periodic 


problem, 


Jy = 2; (6.8) 


where J (without the ig subscripts) is defined by 6.1. Let the matrix A 
be the result of setting by = 0 in the matrix J. Then we rewrite 6.8 in the 
following form: 


Ay = | . (6.9) 


J-1 


Let y 6 be an initial guess. Then we will generate yu from 


(6.10) 


Since A is tridiagonal, we can use Choleski's algorithm to generate yt), 
Thus we form : 


b. 
es ee 
and 
0) (1) 
f, + bry fe Da Oye 
wit) _ 2 che ee ee ag Be Oe yesya Si Ay 
with (6.11) 


1 


(0) () 
(jy... OF as 
J =a a by_yVy-1 


1 
With the sequence {os } thus generated on the forward sweep, we generate 


1 
the sequence yf : on the backward sweep: 


4} 
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A(1) (1) 


= O7 ; 


A(1) (1) A(1), 
Vay Ope Ue s 


(6,12) 
(1) 


A 
tT Y2Y3 3 


(1) 


AN 
tT ViV2 


mW) 2 lt 


Q) _ Q) 
1 


1 
At this point, we shall modify the definition of the sequence {9\ . We de- 


1 
fine the sequence fy ‘ by 


1 1 
yy" = 9° for j a Bs Os awe 


and 


1 f, by 1 1) by f(a 0 
0 eB sag «ale ny AGG 8) 
(2) 


Thus the modification consists of requiring y,; to satisfy the first equation 
of the system 6.9 with the latest estimate available for yy. 


With {,,()} thus generated, the process is repeated, thereby gener- 


ating wo baen ge) The iteration is complete when the last equation in the 


? 


system 6.9 is satisfied to a sufficient degree of accuracy. Hence we ter- 
minate the iteration when 


(2) (£) (2) 


|-byy, ran by_yVj-1 a a a fy (6.13) 


<. e . 
ery 


The code actually stores {+} rather than {e.}; however, the code does 


form Yy, defined by 


] 
Y7 = : (6.14) 
: =s by _Vy-1 


From this, e; can be found and used in Eq. 6.13. Note that this error cri- 
terion is on one component of the residual and thus has nothing to do with 
the rate of convergence of the iterative process. 
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The periodic iterations differ fundamentally from up-scattering iter- 
ations in the following sense. We saw in Section IV that for a source calcu- 
lation (and the same thing is true for a kerf calculation) it does not matter 
how many up-scattering iterations are performed per sweep through the 
channels. The iterative process will still converge to the same solution 
although the rate of convergence will be affected (see Appendix A, Sample 
Problems 8 and 9). The situation with regard to periodic iterations is 
entirely different. Within each channel and for each group, we must find the 
solution, or at least an approximation to the solution, of the system 6,2, 
However, this solution is the solution to a periodic system. The iterates 
y 1), y (2), § dates y (4) we generate to approximate the periodic solution are them- 
selves solutions of a nonperiodic system. Moreover, the higher levels of 
iterations (up-scattering or outer) will not modify the lack of periodicity. 
Thus the periodicity of the solution to the source problem or kerf problem 
is determined by the systems 6.2. From these we must obtain essentially 
periodic solutions (see Appendix A, Sample Problems 6 and 7). 
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VII. INPUT DATA PREPARATION 


In general, ANL-CANDID conforms to the Argonne Standard Reactor 
Code input as specified in ANL-7194.' A type 100000 identification card is 
followed by a set of CANDID code-dependent cards (Type }). 


A. Type 1 Cards 


These cards contain data peculiar to ANL-CANDID: user options, 
problem constants, convergence criteria, criticality search specifications, 
and buckling and volume source data, All data values listed below take the 
default options noted, unless specifically input otherwise. The card type is 
punched beginning in column 1, CANDID reads the cards into one of three 
arrays (entitled INTERFCE, IOPTIONS, and B) depending on column 2, 


1. Type 11 cards for user Options. A blank means the option will 
not be exercised; any nonblank character causes execution of the option. 
Standard problems have nonblank columns 8 and 10. 


INTERFCE 

Subscript Columns Contents 
1-2 Must contain 11. ; 
3-6 Undefined. 

] 7 Process all input cards, but do not execute 
the problem (input check only). 

2 8 Print all fluxes (normally selected). 

9 On-line output (for programmer use only). 

4 10 Print interface fluxes and currents (normally 
selected). 

5 11 Execute slower critical search procedure 
(i.e., run each intermediate k calculation to 
full, specified convergence). This option is 
normally blank, implying the program will 
decide when a control change is to be made. 

6 12 Do not use Chebyshev acceleration during a 
k calculation. This option is normally blank, 
implying Chebyshev acceleration is to be used 
to allow acceleration to speed convergence. 

7 13 Print modified cross sections with every con- 


trol change during a composition search. 
Normally selected if a composition search. 


8 14 Must be blank. 


INTERFCE 
Subscript Columns 
9 15 
10 16 
1] 17 
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Contents — 


Fission fraction (X) vector is the same for 
all materials in cross-section set and does 
not need homogenization into a (space- 
consuming) matrix for every composition, 
normally selected, if possible,to conserve 
storage. 


Perturbation problem is to be executed. 


Two-tape perturbation (real flux on logical 
unit 49, adjoint flux and original cross sections 
on logical unit 48). Otherwise, the pertur- 
bation program will expect one flux tape con- 
sisting of real and adjoint flux and original 
cross sections. 


2. Type 12 cards for Integer data. Integers must be right adjusted 
and fully contained within their six-column fields. 


IOP TIONS 

Subscript Columns 

121 card 1-3 

4-6 

] 7-12 
2 13-18 
3 19-24 
4 25-30 
5 31-36 
6 37-42 
7 43-48 


Contents 


Must contain 121. 

Undefined. 

Number of energy groups in cross-section set. 
Number of energy groups in problem. 

Number of regions in problem. 


Maximum number of outer iterations desired. 
Normally set to 10000 since program will 
stop on time limit. 


Maximum number of up-scattering iterations 
(for up-scattering problems only). Enter 1 if 
not up-scattering problem. t 


Maximum number of fixed fission source 
iterations (generally not used). 


Maximum number of extrapolations on Keff 
during a search problem (4 assumed). 


Tr up-scattering iterations are performed, the convergence rate of the problem Is increased. However, ex- 
perience has indicated that the increase is not sufficient to offset the cost in time for performing the iter- 
ations, Thus, it is recommended that up-scattering problems be tried without using up~-scattering iterations 
first, and these iterations only resorted to if unusual convergence difficulties arise. 
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IOP TIONS 
Subscript Columns Contents 

8 49-54 Maximum number of interpolations on keff 
during a search problem (20 assumed). 

9 55-60 Print frequency during a search (i.e., number 
of control changes in criticality search be- 
fore indicative printing is performed). Nor- 
mally enter a 1 to imply printing after each 
control change. 

10 61-66 Number of general interfaces (feature not 
available). 
11 67-72 = 0 if new problem. 
= 49 if a restart problem. In addition, input 
final k from previous run on 500000 card. 
122 card 1-3 Must contain 122. 
4-6 Undefined. © 
12 7-12 X% direction mesh-refinement factor (relative 
to mesh of input flux guess or restart flux) 
(1 assumed). 
13 13-18 Y direction mesh-refinement factor 
(1 assumed). 
14 19-24 Perturbation-point output key: 61 for printer; 
62 for card punch; 0 if point perturbation out- 
put not desired; another value (1-49) if out- 
put is desired on user logical unit (0 assumed). 
15 25-30 Maximum number of periodic iterations when 
doing an r@ full-circle problem. Normally 
set to 2. 
3. Type 13 cards for Decimal data. All numbers must have decimal 
points. 
B Subscript Columns Contents 
131 card 1-3 Must contain 131. 
3-6 Undefined. 

] 7-18, Desired final value of kes¢ when doing a 
critical search (1.0 assumed). 

2 19-30 Upper bound on search-control parameter 


(X) (0.0 assumed). 


B Subscript Columns Contents 


3 31-42 Second guess (X,) of search-control param- 
eter (0.0 assumed). (First guess on 
Type 500000 card.) 


4 43-54 Lower bound on search-control parameter 
(0.0 assumed). 


5 55-66 Estimate of dk/dx (used to calculate X; if 
nonzero). Normally used, if k(X) is approxi- 
mately known, to eliminate the second control 
pass required to compute dk /dx. Estimate 
may also be available from previous restart 


runs. 
132 card 1-3 Must contain 132. 
3-6 Undefined. 
6 7-18 Flux-convergence criterion on outer bounds 


(k, k) of k (0.001 assumed). 
7 19-30 Not used, 


8 31-42 Smallest flux used in computing monitor in- 
formation (107° assumed). 


9 43-54 Agreement required in neutron balance 
equation (0.001 assumed). 


10 55-66 Estimate of initial kes in critical search 
problem (1.0 assumed). 


133 card 1-3 Must contain 133. 
3-6 Undefined. 

1] 7-18 Extrapolation parameter for fixed fission 
source iteration (1.0 assumed). Normally 
not used, 

12 19-30 Up-scatter iteration convergence criterion 


(0.001 assumed). Normally not used. (See 
footnote on p. 45.) 


13 3]-42 Flux convergence criterion on sum of flux 
differences! (0.001 assumed). 


For a k-calculation, three convergence criteria are required on the outer iterations: 
1. Outer bounds of k: k < k'") < k (132 card) 


1/2 
a 
2, Sum of flux differences: > On bic i) (133 card) 
cig Pig 


3. Ker differences: k!™ - x(N-1) (Type 500000 card) 


We recommend that these all have the same order of magnitude, ranging from 10-4 to 10-7, with a "normal" value of 10-5. That is, a convergence of 1074 may be desired 
for some purposes, but does not generally yield sufficient accuracy. A convergence of 107 yields sufficient accuracy, but at a large cost in computer time. Hence, this 
convergence is not generally used untess the results are to be input to the perturbation code where small perturbations are to be made. 
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48 


B Subscript 


14 
15 


16 


17 


18 


4, 


Columns 


43-54 
55-66 
1-3 
3-6 

7-18 


19330 


31-36 


Contents 


Not used. 

Asymptotic inverse reactor time period (a). 
Must contain 134. 

Undefined, 


Initial estimate of outer iteration convergence 
rate used in computing Chebyshev extrapola- 
tion coefficients (0.995 assumed). Normally 
set to 0.995. (See page 37.) 


Periodic iterations convergence criterion 
(rQ full circle). Normally set to 0.0001. 


Interval shift factor for Chebyshev accelera- 
tion (0.0 is assumed, implying no shift). 
Normally set to 0.0. 


Type 14 and 15 cards for Buckling and Volume Source (decimal) 


data. Data can be input as one constant, as varying with energy only or with 
region only, or as varying both with energy and region number. 


Columns 


1-2 
3-6 


7-18 
19-30 
31-42 
43-54 
55-66 


67-69 


noted below. 


B. Type 2XXXXX Cards 


Contents 


Must contain 14 for Buckling data and 15 for Volume Source 


data. 


Region number, given only when data are to be distributed 


into regions. 


Buckling or Volume Source data, as named in columns 1-2, 
given on as many cards as necessary to cover all energy 
groups for this region, before beginning next region. 


Energy-group number of first group on this card, given only 
when data are to be distributed by energy groups. 


The remaining cards in the input deck are standard, except as 


The Type 210000 cards for geometry specification must be consistent 
with mesh of the input flux guess. If a problem is being restarted using an 


old flux tape, the number of mesh points and the region structure given on the 
Type 2 cards must correspond to the tape structure. Mesh "multiplying" in 


X or Y directions is accomplished by setting the values of the X and Y mesh- 
refinement factors (IOPTIONS (12) and (13) on the 122 card). The factor in- 
put is the number of mesh points to be made from each input mesh point. The 
problem is then run using the new, "full" mesh. Subsequent restart problems 
to be run using the new mesh will require a new set of Type 2 cards con- 
sistent with this mesh (and the removal of mesh-multiplying factors 
IOPTIONS (12) and (13)). That is, mesh refinement is used only once for 
each desired mesh size without redefining the mesh in the Type 2 cards. 


The Type 220000 cards for the ratio-method geometry specification 
is not implemented. 


C. Type 3XXXXX Cards 


Column 60 is ignored; i.e., general boundary conditions are not 
implemented. 


On the remaining Type 3XXXXX cards the only restriction on card 
ordering is that the sides be input one at a time, until all regions and groups 
on a side are covered, with any ordering of data within a side allowed. 


D. Type 4XXXXX Cards 


Both Type 400000 and 410000 cards must always be present. That 
is, direct isotope-to-composition homogenization is not implemented. 


E. Type 5XXXXX Cards 


On the Type 500000 card, the value for 6 (columns 13-24) is used as 
the convergence criterion on differences between k values; for generation 
n, the criterion is satisfied if jk(n) - k(n-1)| < 6 (see footnote on p. 47). The 
initial guess supplied in columns 46-57 is used as an initial value for k in 
a flux calculation, or for uw in a source calculation (blank or 0 signals the 
default option of 1.0). The k corresponding to a restart flux must be used 
for a restart problem. If the problem is a criticality search, the initial 
guess, i.e., the control parameter (Xj), is placed in this field, and the 
initial k is put on the 132 card. The remaining Type 5XXXXX cards are 
implemented as specified. All cards with region numbers must have them 
in increasing order. Energy-dependent buckling modifiers can be input 
using the concentration change format, with columns 2-6 containing a region 
number, if desired, and columns 19-24, 37-42, and 55-60 containing group 
numbers corresponding to modifiers in columns 7-18, 25-36, and 43-54. 


F. Type 6XXXxXxX Cards 
On the Type 600000 card CANDID requires the initial guess to bea 


flux guess andnota source; therefore column 7 must benonblank. Spectrum 
averaging is not performed. External source mustbe read from Type l5 cards. 
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G. Type 7XXXXX Cards 


These cards are as specified in ANL-7194. In addition to the normal 
use of cross-section set modification,! these cards are used to input cross 
sections for use in perturbation calculations.® 


TANL-CANDID expects microscopic cross-section data on a magnetic tape prepared by the XLIBIT Code. 
An XLIBIT cross-section set must be available on tape for each problem even if all cross-section data are 


to be input via the 700000 cards, 
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VIII. EXECUTION CARD DECK ARRANGEMENTS 


| A. Normal Execution 


The control cards and deck structure necessary to execute a normal 
ANL-CANDID run are (§ punches are assumed in column 1): 


JOB, User Identification 

MOUNT XLIBIT Cross-section Tape on LUN 44 

MOUNT CANDID AOV Tape on LUN 47 (AOV means absolute overlay) 
SCRATCH TAPES on LUN 42, 43, 45, 55 (SYSTEM scratch unit) 
EQUIP,44=(XLIBIT Cross-section tape name, Version) RO, SV 
EQUIP,3=45 (Equivalent reference saves one scratch tape) 
EQUIP,47=(CANDID AOV tape name, Version), RO, SV 

LOAD MAIN, 47, Time estimate, Print estimate 

Data 

END OF PROBLEM (column 3) 


B. Program Compilation 


The volume of source program cards (some 20000) in ANL-CANDID 
prohibits the handling of source decks directly. Thus, the cards are put on 
tape along with SCOPE cards denoting the end of overlays. The following 
SCOPE control cards (in addition to special mounting instructions) will com- 
pile the source code from LUN 10 (unlabeled) and prepare the binary load 
and go (LGO) tape on LUN 11. Note that CANDID is compiled as a one-bank 
code, 


MOUNT TAPE (ANL-CANDID source tape reference number) ON 
LUN 10 

EQUIP,10=SV ,** 

EQUIP,11=(CANDID LGO tape name, Version), SV 

FILE,11 

MAIN ,47 

FILE END Main Control Section 

FIN, 4=11,1=10;R.*: 

COMPASS,L,X=11,1=10 

FILE,I1 

OVERLAY ,47,1 

FILE END 

FIN ,g,.5=11,1=10,R,%. 

FILE,11 

OVERLAY ,47,2 

FILE END 

FTN,L,X=11,1=10,R,*. 


Overlay 1: Processes card Types 1 and 2, 


Overlay 2: Processes card Types 3-5. 
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FILE,11 

OVERLAY ,47,3 

FILE END 
FUN, x21 11=10.8;*. 
FILE,1I] 

OVERLAY ,47,4 

FILE END 
FTN,L,X=11,I=10,R,*. 
FILE,11 

OVERLAY ,47,5 

FILE END 
PIN: % == 10,R.*. 
FILE,!1 

OVERLAY ,47,6 

FILE END 
PUN xsi ia10-R 
FILE,11 

OVERLAY ,47,7 

FILE END 

ETN =) i= 10.R st: 


Overlay 3: Creates pseudo cross-section 
library tape from XLIBIT tape. 


Overlay 4: Processes card Type 7 if 
they exist. 


Overlay 5: Performs search and k 
calculations. 


Overlay 6: Output edit. 


Ks ee ee 


Overlay 7: Perturbation. 


C. Program Modification Using LGOEDIT 


Additions, corrections, or deletions to the load-and-go tape are made 
using the program LGOEDIT (see Appendix B). The card deck structure 
listed below will accomplish the following (assuming the existence of an old 
load-and-go tape): 


1, Compile the indicated source programs and place them, together 
with the associated LGOEDIT control cards (REP,DEL,INS), on the merge 
tape (LUN 45). 


2. Select the first program name (NAME1) from the control card 
on the merge tape, and search the old LGO tape (LUN13) until NAME! is found. 
Preceding subroutines are simply copied from LUN1I3 onto the new LGO 
tape (LUN11). The instruction on the control card is executed as described 
in Appendix B, and LGOEDIT proceeds to the next subroutine on the merge 
tape. Note that LGOEDIT does not space the old LGO tape, indicating that 
the ordering on the old LGO tape must be preserved on the merge tape. 


MOUNT TAPE (Old LGO tape) on LUN13 
EQUIP,13 = (Old LGO tape name, Version), RO,SV 
EQUIP,11 = (New LGO tape name, Version), SV 
FILE, 45 

REP NAME 1 


FILE END 
FTN,L,X=45,*. 
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Source for Subroutine NAME] 
SCOPE 


FILE,45 

DEL NAME2 
FILE END 
FILE,45 

INS NAME3 
FILE END 
FTN,L,X=45,*. 


Source for Subroutine NAME3 
SCOPE 


Binary Deck for LGOEDIT 


RUN 
11 | 13 | 45 (Tape numbers as input to LGOEDIT columns 1-6) 


D. Modification Followed by Execution 


The normal execution structure (A) may be merged with the program 
modification procedure (C) to produce a run in which modification is followed 
by execution. The structure is as follows: 


Normal control cards (A) 
Modification structure (C) 


LOAD, 11 (Prepares AOV Tape on LUN 47) 


RUN 
Problem Data 


E. Restart 

A restart feature is provided for the following reasons: 

1. Machine failure usually in the form of tape parity error. 

2. Time estimate exceeded. 

CANDID checks internally for an error condition, and when one is 
encountered, the current flux vector is written, as indicated in Fig. 4, onto 
a flux save tape (LUN 48). Normal termination will also result in the flux 
being written on LUN 48. Thus, to save the flux, the user must insert the 


following control card: 


EQUIP,48=(Flux save tape name, Version), SV. 
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Original Composition 
Cross Sections 
Real Flux Adjoint Flux | (Used by Perturbation) 


Cross sections 
by composition 
and group 


Same format 
as Real Fiux 
(Dimensions 


First channe! 


may differ type cross 
. during mesh sections 
First refinement) 
Channel 
Flux 
Lth channel! 
type cross 
sections 
Second 
Channel 
Flux 
Ith 
Channel 
Flux 


Fig. 4, Restart Tape Format 


To restart the problem from this point, the user must remount the 
tape on the restart unit (LUN49) by inserting the following SCOPE cards: 


Mount tape (Old Flux tape Identification) on LUN 49 
EQUIP,49=(Old Flux Tape Name, Version), RO,SV 


The original input-data deck is used with the exceptions (1) the 
number 49 must be punched on the 12 card (see Section VII), and (2) the 


55 


mer where n is the last iteration of the previous run, must be punched 
on the 500000 card (or on the 132 card if a search problem). Normal exe- 
cution will then proceed with the flux guess taken from LUN 49 instead of 
the standard input unit. The new flux may again be saved on LUN 48. (See 
Section IX.) 


F. Back-to-back Runs 


The stacking of runs back-to-back is provided, but is not usually 
done except when perturbation problems are run.® The reason for this is 
that two-dimensional problems of any size require so much machine time 
that they are run in stages. (See Section IX.) 


A complete exit to the SCOPE System is made following each End of 
Problem card. Thus, the program must be reloaded for each problem. This 
is done to allow for any available mounting or dismounting type instructions 
to be recognized between problems. Thus, in stacking problems, the deck 
structure is: 


Normal execution control cards (A) 

Program modification structure, if any (C) 

LOAD, RUN 

OR 

LOADMAIN ,47 

Data Problem 1 

END OF PROBLEM 

Tape-handling instructions (i.e., unloading old flux tape, mounting 
new flux tapes for next problem, etc.) 

LOADMAIN ,47 

Data Problem 2 

END OF PROBLEM 


Tape-handling instructions 
LOADMAIN ,47 

Data Problem n 

END OF PROBLEM 

End of File 


G. Sample Deck Structure 


A sample deck structure including program modification, restart, 
and back-to-back problems might appear as follows: 


JOB, User Identification 
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MOUNT OLD LGO Tape ON LUNI1 

MOUNT XLIBIT Cross-section Tape on LUN 44 
MOUNT Old Flux Tape for Problem 11 on LUN 49 
SCRATCH Tapes on LUN's 42,43,45,55 
EQUIP,11=(NEW CANDID LGO,L),SV 
EQUIP,13=(OLD CANDID LGO,),SV 

EQUIP,3=45 

EQUIP ,44=(XLIBIT Cross-section Tape Name,Version),SV 
EQUIP,47=(CANDID AOV,2),SV 
EQUIP,48=(NEWFLUX1,1),SV 
EQUIP,49=(OLDFLUX1,1),SV 

FILE ,45 

REP NAME! 

FILE END 

FTN,L,X=45,*. 


Source deck for Subroutine NAME] 
SCOPE 
Binary deck for LGOEDIT 


RUN 

131145 (LGOEDIT input card) 

End of File 

Release, 13 (Special ANL Control Card which unloads 13) 
Load, 11 

Release, 11 

RUN 


CANDID Data Deck for Problem 1 


End of Problem 

End of File 

Release, 48 

Release, 49 

Mount old flux tape for Problem 2 on LUN 49 
EQUIP,48=(NEWFLUX2,1) ,SV 
EQUIP,49=(OLDFLUX2,1) ,SV 
LOADMAIN ,47 


CANDID Data Deck for Problem 2 


End of Problem 
End of File 
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IX, NORMAL EXECUTION PROCEDURES 


A two-dimensional diffusion problem is usually run in stages (a 
sequence of computer runs) for two reasons: 


1, The probability of error in problem input data (even if it gets 
by the input processor, which checks for consistency only) is significant. 
Thus, we would like to be assured that the problem is proceeding in the 
"right" direction before spending a significant amount of computer time. 


2. Even when we are assured that the problem is progressing 
smoothly, it is worthwhile to continue in stages of less than 60 min of com- 
puting time to measure convergence rates, progress, etc. 


The stages mentioned above are of two types: 


1. Mesh size: i.e., starting with a coarse mesh, running the 
coarse-mesh problem, and refining the mesh using the coarse-mesh flux 
as an input flux guess for the refined problem. 


2. Computer time: i.e., putting a time limit on execution of 
something less than 60 min, coming off the computer to measure progress, 
and restarting and continuing the problem from that point. 


Experience has shown that a 20-group, 40 x 40 mesh may take 3 to 
4 hr of computing time (depending on accuracy requirements and the rate of 
convergence of the problem). Hence, it is wise not to be hasty in expending 
this amount of time only to find (possibly because of an error or incorrect 
assumptions) that the results are not what was expected. The problems are 
generally run as follows: 


1, The mesh is reduced or "coarsened" to a problem that will run 
in 20 or 30 min (e.g., 10 x 10 x number of groups); the results are expected 
for "ballpark" accuracy. Experience indicates that this result is certainly 
within 10% of the final answer and may be as goodas 1%, Errors in con- 
figuration should show up at this time. At the end of the job, the flux tape 
is saved as indicated under Execution Card Deck Arrangements, 

Section VIII. 


2. After step 1 is completed and the refined problem is to be at- 
tempted, the mesh-refinement feature is used. This amounts to placing a 
multiplication factor (integer) for the X and/or Y direction on the 122 card 
(see Section VII). For example, a 2 placed in the field for the X direction 
will double the mesh along this axis. (The coarse problem must be planned 
so that an integral multiple of the coarse mesh will give the final desired 
mesh.) CANDID now automatically spreads the coarse mesh to adjacent 
points in the fine mesh, as indicated in Figs. 5 and 6. The resulting flux 
guess is then used to start the problem. Refinement can take place as often 
as desired along either axis. 
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| 1 ‘ ‘ i] ] ] | 
r, to 3 4 5 Ve v7 Veg 


Coarse mesh: O Flux point 
Refined mesh: X Flux point 


The orginal 3 x 4 mesh (points O) 6 x 8 (points X) 


Fig. 5. Mesh Refinement* (Doubling in r and z directions) 


r 
ro ry lo rs V4 
| t i) i] 4 1 i] $ 
r, Vo 3 a rs Ve V7 Teg 
Note: An interpolation scheme would give a better approximation, 
but one has not been incorporated into ANL—CANDID. 


Fig. 6. Typical Flux Shape for Starting the Refined Problem* 


**Unprimed mesh lines denote coarse mesh. 
Primed mesh lines denote refined mesh. 
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X. COMPUTING REQUIREMENTS AND EXECUTION TIME 
A, Machine: CDC 3600 with one-bank 32K fast memory. 
B. Operating System: SCOPE 


C. Tape Units: Normal Execution: Four scratch tapes, cross~section tape, 
AQV tape. 


Program Modification: Four scratch tapes, cross-section 
tape, old LGO tape, new LGO tape, AOV tape. 


Restart: Four scratch tapes, cross-section tape, AOV tape, 
old flux tape, new flux tape. 


D. Execution Time: Execution time depends on the convergence criterion 
| and the problem convergence rate. Typically, a 
40 x 40 mesh, 20-group problem will run 3-4 hr. A 
10 x 10 mesh, 20-group problem will run 20-30 min. 
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XI, RESTRICTIONS AND LIMITATIONS 


1. Roughly, the largest problem that can be handled is 20 groups 
by 40-45 points in a channel (i.e., Y or Z or @ direction). Any number of 
points can be handled in the R direction if one has enough money. 


2. The number of compositions cannot exceed the number of ma- 
terials in the cross-section library. "Dummy" isotopes should be inserted 
in the cross-section set if this situation arises. 


3, Composition cross sections cannot be directly homogenized 
from isotopes. The material level cross sections must be inserted even 
though they may be decimals (i.e., volume fractions are 1.0). 


4, A source guess cannot be input. Only a flux guess is allowed. 
5. The distributed source for a source calculation must be input 


by region and group on the Type 15 cards. CANDID will not accept an ex- 
ternal source input via the 610000 cards. 


6. Generally, the output is erroneous for adjoint calculations, but 
the flux printout is correct. 


XII. SUMMARY OF REQUIREMENTS 


Execution of ANL-CANDID requires familiarity with the following 
documents, in addition to ANL-7305: 


1. Code independent input and cross-section library tape prepara- 
tion are described in the following documents: 


M. Butler and H. Greenspan, General Input Specifications for 
ANL Reactor Programs, ANL-7194 (April 1966). 


S. D. Sparck, XLIBIT: An ANL Cross-section Library Code, 
ANL-7112 (Feb 1966). 


toe To make modifications and/or additions to the CANDID load- 
and-go tape, the use of LGOEDIT is helpful: 


LGOEDIT2, reproduced as Appendix B by permission of the 
authors. 


3. For perturbation calculations, the following publication is 
required: 


G. K. Leaf and A. S. Kennedy, PERC, A Two-dimensional 


Perturbation Code Based on Diffusion Theory, ANL-7304 
(May 1967). 


4, Finally, if available, the Control Data Corporation documen- 
tation of CANDID2D may be helpful: 


CANDID2D, A Two-dimensional Neutron Diffusion Program, 
Control Data Corporation (1966). 


The following items are necessary to run ANL-CANDID: 


1. Source (BCD FORTRAN Card images) tape. 

Load-and-go tape prepared by compiling source tape. 
Absolute overlay tape prepared from the load-and-go tape. 
XLIBIT generated cross-section library tape. 

Proper SCOPE control cards. 

Input deck which conforms to ANL-7305 and ANL-7194. 


NY Oo WwW Ff W W 


Lots of luck!! 
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APPENDIX A 


Sample Problems 
The sample problems in this appendix reinforce the conclusions 
drawn in the body of the report. Hence only relevant parts of each output 


listing are included. For completeness, sample problem 1 is listed in 
its entirety. 


For reference, the tests for convergence are listed below: 


2 1/2 
1. Sum of flux differences = \ Oe - gin} , 
"3 1g 1jg 
1,J,8 
where 
(i,j) are the point indices, 
g is the group index, 
and 
n is the outer iteration index. 
2. Bounds of kKogp = k- &:, 
where 
(n) 
a Pig 
k = max 
ijg gard) 
jg 
and 
gin) 
1jg 
k= mins 
— “ijg (n-1) 
“13g 


(n) | (n-1) 
keep - Kege [> 


S. Kare difference = 


(n) 
kerf - Kgesired 


4. Criticality search = 


yy go) pee) 
; JS JE 
5. Up-scattering = JB 
sy (1) (0) 
LlPig ~ Pig 
J& 
6. Periodic = -byy" - by yy) + esy\”) 7 


1. Sample Problem 1 
a. Description 
(1) Problem Type 
Real kere calculation without Chebyshev acceleration. 
(2) Configuration 
(a) Geometry 
rz 
(b) Region Definition 


The reactor consists of three regions as follows: 


Region No. 1: Core region composed of uranium 

2 238 : ; 
(U 35 and U” ), stainless steel, molybdenum, zirco- 
nium, iron, nickel, chromium, and sodium. 


Region No. 2: Control- and safety-rod region com- 
posed of uranium (U*? and U”*®), molybdenum, nio- 
bium, iron, nickel, chromium, and sodium. 


Region No. 3: Radial-blanket region composed of 
uranium (U2 and U8), plutonium (Pu2’), iron, 
nickel, chromium, and sodium. 


(c) Mesh Definition 


r direction: 27 points 


z direction: 20 points 


63 


64 


(d) Boundary Conditions 
Left: ¢! = 0. 
Right: ¢ = 0. 
Bottom: @¢! = Q. 
Top: @ = 0. 
(e) Number of Energy Groups: 2 


(3) Convergence Criteria 


= = 5 
keg difference = 10 ”, 


ke¢p bounds = 107°. 


Sum of flux difference = 107°. 
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2. Sample Problem 2 
a. Description 
(1) Problem Type 
Real kg¢¢ with Chebyshev acceleration 
(2) Configuration 
Same as Sample Problem 1 


(3) Convergence Criteria 


Same as Sample Problem 1 
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3. Sample Problem 3 
a. Description 
(1) Problem Type 


Dimension search varying the z direction uniformly to 
obtain Kage =-'102 5». 


(2) Configuration 
Same as Sample Problem l. 
(3) Convergence Criteria 


kere difference = 107°, 
Kare bounds = 10-5, 
Sum of flux differences = 107°. 


kage = 1,025) = 5% 10°°, 
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4, Sample Problem 4 


a. 


(1) 


(2) 


Description 


Problem Type 


Source calculation with Chebyshev acceleration and 


= 1/0.9. 


(a) 


(c) 


(d) 


Configuration 


Geometry 


r@ (quarter circle) 


Region Definition 


The reactor consists of 18 regions as follows: 
Region No. 1: Internal thermal column (water). 
Region No. 2: Stainless steel shroud. 


Regions No, 3-11: Fuel regions composed of enriched 
uranium, stainless steel, and water. 


Region No. 12: Control-rod region composed of zirco- 


nium and water. 


Regions No. 13-15: Reflector composed of beryllium, 
aluminum, and water. 


Region No. 16: Outer reflector composed of aluminum 
and water. 


Region No. 17: Iron vessel wall. 


Region No, 18: Water. 
Mesh Definition 


r direction: 33 points. 
Q@ direction: three points. 


Buckling (B*) = 0.052360. 


Boundary Conditions 


Left: g' = 0. 
Right: @ = 0. 
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Bottom: ¢! = 0. 


Top: @' = QO. 
(e) Source 
Neutron source density in region 15: 1.0. 
(f{) Number of Energy Groups: 16. 
(3) Convergence Criterion 


Sum of flux difference = 107°. 
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5. Sample Problem 5 
a. Description 


(1) Problem Type 


Source calculation without Chebyshev acceleration and 


i= 1/0.9. 
(2) Configuration 
Same as Sample Problem 4. 


(3) Convergence Criteria 


Same as Sample Problem 4. 
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6. Sample Problem 6 


a. 


Description 


(1) Problem Type 


(2) 


(3) 


Real Kore calculation 


(a) 


(b) 


(c) 


(d) 


(e) 


Configuration 


Geometry 


r@ (full periodic) 
Region Definition 


The reactor consists of a core region composed of 
: 23 : ; 
uranium (U > and U): iron, and aluminum. 


Mesh Definition 


r direction: 20 points 


0 direction: 20 points 


Boundary Conditions 


ett: ol = 0; 
Right: © = 0, 
Bottom: periodic. 


Top: periodic. 


Number of Energy Groups: 2. 


Convergence Criteria 


Kate difference = 107°. 


Koff bounds = 107°. 


Sum of flux difference = 107°, 


Periodic = 107%, 


Maximum number of periodic iterations = 25. 
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7, Sample Problem 7 


a. Description 


(1) 


Problem Type 


Real Kage calculation 


Configuration 


Same as Sample Problem 6 


Convergence Criteria 


kage difference = 107°, 

Kate bounds = 107%. 

Sum of flux difference = 107%. 
Periodic = 107%. 


Maximum number of periodic iterations 
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8. Sample Problem 8 
a. Description 
(1) Problem Type 
Real kg¢r calculation 
(2) Configuration 
(a) Geometry 


xy 
(b) Region Definition 
The reactor consists of three regions as follows: 


Region No. 1: Internal thermal column (water). 


Region No. 1: Core region composed of enriched 
uranium, stainless steel, and water. 


Region No. 3: Reflector composed of beryllium and 
water. 


(c) Mesh Definition 


x direction: nine points 
y direction: nine points 


Buckling (B*) = 0.0027416 
(a) Boundary Conditions 

Left: ¢' = 0. 

Right: @ = QO. 

Bottom: @¢ = 0. 

Top: @ = 0. 


(e) Number of Energy Groups: 18. 
(3) Convergence Criteria 


ket difference = 107°. 

keff bounds = 107°, 

Sum of flux difference = 107°. 
Up-scattering = 107°. 


Maximum number of up-scattering iterations = Q. 
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9. Sample Problem 9 
a. Description 
(1) Problem Type 
Real kore calculation 
(2) Configuration 
Same as Sample Problem 8 
(3) Convergence Criteria 


kare difference = 10°? 
Kore bounds = 10: =. 
Sum of flux difference = 107°. 


Up-scattering = 107°, 


Maximum number of up-scattering iterations 
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APPENDIX B 


LGOEDIT2 
1, Identification 
a. Title. Load-and-go tape edit 
b. Language. FORTRAN-63, one COMPASS subroutine 
c. Machine. 3600 
d, Programming Sanford Elkin,t D. B. Taylor, A. L. Rago, 
and Design. J. Zapatka, G. J. Duffy 


2, Purpose 


LGOEDIT2 provides a facility for editing and updating load-and-go 
tapes for the 3600 by means of deleting, replacing, and inserting binary 
decks and control cards at specified positions on the tape. LGOEDIT2 is 
expected to be useful primarily when the LGO tape contains MAIN, 
OVERLAY, and SEGMENT control cards, or when the order of loading 
is important. 


3. Usage 


LGOEDIT2 creates a new LGO tape by editing an old LGO tape with 
information from an edit tape (the edit tape can be the standard input unit). 
LGOEDIT2 first reads a card from standard input using the format 31¢. 
Colums 1 and 2 must contain the logical unit number of the old LGO tape, 
colums 3 and 4 the logical unit number of the new LGO tape, and columns 
5 and 6 the logical unit number of the edit tape. (All must be numbers 
from 1 to 49, except that columns 5 and 6 can contain 60.) All further 
information is obtained from an edit tape. 


LGOEDIT2 recognizes five control cards. These are: 


Column 1-4 5-8 9-16 
§DEL NAME 
g§REP NAME 
SINS | NAME 
SBYP | ; NAME 
{Dic NAME 


T Control Data Corporation. 


de 


where NAME is the BCD name of the subroutine of the LGO tape, and it 
must be left justified beginning in column 9. Columns 5-8 are not used. 


LGOEDIT2 reads a control card from the edit tape and copies the 
old LGO tape onto the new one until it finds subroutine NAME. It then does 
the operation required by the control card and returns to the edit tape to 
read another control card. When it finds an end of file on the edit tape, 
editing is completed, and the remainder of the old LGO tape is copied 
onto the new one. 


The control cards cause LGOEDIT2 to perform the following 
actions: 


SDEL will delete subroutine NAME. 


§REP will replace subroutine NAME with binary decks from the 
edit tape, until another control card is read. Note that NAME 
may be replaced by many subroutines. 


gINS will insert binary decks from the edit tape onto the new LGO 
tape after subroutine NAME has been copied, Note that many 
decks may be inserted after NAME. 


§BYP will copy the old LGO tape onto the new LGO tape, up to and 
including NAME. The reason for this is that on an overlay tape, 
the same subroutine name may occur in several overlays or 
segments. If we desire to edit some of these subroutines, sBYP 
allows us to bypass an earlier occurrence of NAME. 


5DLC will delete all cards between the IDC card of subprogram 
NAME and the final card of the subprogram preceding NAME. 


This DLC card was installed mainly to remove or replace loader 
control cards (e.g., bank, main, overlay, segment) or octal 
corrector cards. 


LGOEDIT2 never looks at the names of the subroutines on the 
edit tape, but only at the NAMES on the control cards. 


4. Restrictions 


The NAMEs referenced on the edit tape must appear in the same 
order as on the old LGO tape, because LGOEDIT2 will not backspace the 
tapes to search for missing names. Both the edit tape and the old LGO 
tape must be terminated by an end of file. We cannot follow a REP card 
with an INS card with the same NAME on both. 
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5, Timing 
A large LGOEDITZ job was timed with the following resulting sta- 
tistics: 
a, LGOEDIT2 portion of job: 
Number of DLC control cards = 2 
Number of DEL control cards = 5 
Number of REP control cards = 8 
Number of binary subroutines transferred = 143, 
LGOEDIT2 time = 1 min 35 sec (22% of job time) 


b. FORTRAN called nine times to compile 11 subroutines (1650 
source cards). 


FTN time = 2 min 50 sec (40% of job time) 
c. Translation of LOAD/GO tape to Overlay tape. 
Translation time = 2 min 44 sec (38% of job time) 


Total time for all three runs in job = 7 min 9 sec 


6. LGOEDIT2 Output Edit 


LGOEDIT2 will produce lists of the subprogram names on old and 
new LGO tapes, along with a record of the edit procedures carried out 
and attempted. 


7. Typical LGOEDIT Job 


?Job,999,204051,25 

Acct. card. 

3/Mount Tapes L5882(44) and L8878(49) 

SE QUIP, 20=60 

EQUIP, 21=61 

SE QUIP, 29=(SNARG:-OVERLAY:29,1),SV 

'— QUIP,39=(SNARG:-OVERLAY:-LOAD-AND.GO,1),SV 


SE QUIP,44=(SNARG-X LIBIT- DEBUG, 6)RO,SV 
'e QUIP,49=(SNARG:-OVERLAY- LOAD: AND.GO),RO,SV 


SFILE,14 
jSREP...--SNARGID 
SFILE-END 
s5F TN,L,A,X=14. 
SNARGID Source Deck 
SCOPE 
jFILE, 14 
SREP: INITILIZ 
§FILE-END 
SF TN,L,A,X=14 
INITILIZ Source Deck 


SCOPE 


LGOEDIT2 Binary Deck 


JRUN,60,10000,7, 
493914 


SEND-OF’FILE card for end of this run 


§LOAD,39 
§RUN,60,10000,7, 
Data Deck 
Blank Card 
send of File card for this job. 
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To place the following on 
LUN 14: 


jREP:---SNARGID 
SNARGID binary deck. 
SREP..--INITILIZ 
INITILIZ binary deck. 
3(EOF) 


Transfers all binary card records from 49 to 39 replacing only 
those belonging to SNARGID and INITILIZ in the process. These 
two are transferred from 14 to 39 in lieu of their old versions. 
After the LGOEDIT2 run, LUN39 was loaded to run the program. 
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